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This paper presents comparisons of seven propagation codes Tor predicting liner attenuation in ducts with 
How. The selected codes span the spectrum of methods available (finite element, parabolic approximation, and 
pseudo-lime domain) and are collectively representative of the state-ol-art in the liner industry. These codes are 
included because they have two-dimensional and three-dimensional versions and can be exported to NASA’s 
Columbia Supercomputer. The basic assumptions, governing differential equations, boundary conditions, and 
numerical methods underlying each code are briefly reviewed and an assessment is performed based on two 
predefined metrics. The Hvo metrics used in the assessment are the accuracy of the predicted attenuation 
and the amount of wall clock time to predict the attenuation. The assessment is performed over a range of 
frequencies, mean tlow rates, and grazing flow liner impedances commonly used in the liner industry. The 
primary conclusions of the study are (1) predicted attenuations are in good agreement for rigid wall ducts, 
(2) the majority of codes compare well to each other and to approximate results From mode theory for sott wall 
ducts, (3) most codes compare well to measured data on a statistical basis, (4) only the finite element codes with 
cubic Her mile polynomials capture extremely large attenuations, and ( 5 ) wall clock time increases by an order 
of magnitude or more are observed for a three-dimensional code relative to the corresponding two-dimensional 
v ersion of the same code. 


I. Introduction 

H igh levels of engine fan noise radiated from modem turbofan engines threaten to curtail much needed growth in 
commercial air transportation systems worldwide, Engine fan noise is acoustic energy generated at the fan or e\it 
guide vanes of a turbofan engine that radiates forward via the nacelle inlet duct and rearward by propagating through 
the aft nacelle duct and then refracting through the ambient shear layer. Figure 1 shows a schematic of a typical aircraft 
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Fig. 1 Schematic of a turbofan engine nacelle with wall lining. 

engine nacelle. Earlier versions of the turbofan engine were dominated by tones associated with the blade passage 
frequency anti ils harmonics. However, in recent engines the fan blade shapes were modified to incorporate lean and 
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sw eep, such that the blade passage frequency is less dominant. Newer engines have also incorporated increased fail 
bypass to increase the thrust to drag ratio mid fuel efficiency 

One of the most effective means for reducing levels of noise radiating from commercial aircraft engines has been 
the installation of acoustic treatment in the walls of the inlet mid aft fan ducts. Such treatments reduce radiated noise 
by attenuating inlet mid aft fan noise as it propagates through the engine nacelle. Suitable aircraft noise prediction 
codes are required to optimally design the acoustic treatment for maximum attenuation within a specified nacelle 
geometry Time domain methods are a topic of considerable research interest; however, the frequency domain method 
remains the primary basis for the prediction of liner attenuation in the aircraft noise industry. Although engine designs 
incorporate three-dimensional (3D) geometries, a literature survey reveals that until the turn of the 21st century, the 
only numerical method for predicting liner attenuations in 3D geometries was the parabolic approximation. 1 Since 
that time, the finite element method 2 and the pseudo-time domain method 3 have also emerged. In addition, the NASA 
Langley Research Center has two finite element codes 4-6 that were recently extended to 3D. 

Each computer code that has been written based on the above three methods has advantages and disadvantages 
depending upon the problem being modeled, assumptions that arc made, mid computer facilities available to code 
users. Scientists desiring to use codes based on the above three methods are often frustrated because (1) information 
available to code users is often fragmented and incomplete, and much of the theory underlying code development {c.g., 
the assumptions made, the governing differential equations, the boundary conditions, and die numerical method) is 
often not clear, (2) comparison of wall clock time statistics between the various codes are not available because codes 
have not been run for the same problem, using the same computer, and on identical grids, (3) exact analytical and 
approximate attenuations from mode theory are av ailable to compare with code predictions, but these comparisons 
are rare, (4) comparison of code predictions with measured data is at best questionable because experiments are 
seldom designed to prov ide the necessary code inputs (eg., the sound source, grazing flow liner impedance, and duct 
termination conditions), and (5) the measured data contain uncertainties in the inpuLs needed for code predictions; 
however, code users have treated these inputs as deterministic. 

The purposes of this paper are to (1) collect the most commonly used 3D duel propagation codes (non- vortical ;md 
vortical) for predicting liner attenuation under a single umbrella, (2) briefly review the basic assumptions, gov erning 
differential equations, boundary conditions, and numerical method underlying code development, mid (3) perform an 
assessment based upon two predefined metrics. The selected codes span the spectrum of methods currently in use 
(parabolic approximation, finite element method, and pseudo-time domain method) and are collectively representative 
of the current state-of-art in the liner industry The two metrics used in this assessment are the accuracy of the predicted 
attenuation and the wall clock time required to predict the attenuation. 

The following section (Le., section II) describes the physical duct acoustic problem of interest and the basis for 
code selection. Sections III and IV present the basic assumptions, the governing differential equations, duct boundary 
conditions, and numerical method underlying code development for both the non- vortical mid vortical ducL propagation 
codes. The method used to predict liner attenuation in the presence of grazing flow is described in section V, Results 
of the code assessment are presented in section VI, and a brief summary along with the primary conclusions of this 
study are presented in the final section of this paper, section VII. 


II. Statement of Problem and Basis for Code Selection 

F igure 1 shows a sketch of a typical 3D duet geometry. The duet carries a non-uniform, subsonic mean flow from 
left to right. A Cartesian coordinate system is used, in which the axial coordinate is z and the two transverse 

coordinates arc x and y. The source plane is denoted by an imaginary plane, V s - Two types of acoustic propagation 
are considered For inlet propagation, the sound propagates upstream of the fan against ihe flow and exits the duct 
as shown in the upper portion of Fig. 2. In this case, the duct termination is located at the entrance to the inlet and 
is denoted by the plane, IV The second ty pe of propagation is aft-duct propagation, in which the sound propagates 
downstream of the fan face and exits the duct via the alt -fan nacelle duct as shown in the lower portion of Fig, 2. For 
aft propagation, the termination plane, V !y is located downstream of the fan face in the exhaust duct. In both inlet 
and aft duct propagation, the computational v olume will be denoted by £2 and the perimeter of the duct containing the 
physical duct walls will be denoted by F u -. Let the surface bounding Q be denoted by F, such that the union of F^, F,, 
and F w is F 

r = r s ur t ur w ( 1 ) 

The perimeter of the duct is lined with acoustic material that is assumed to be locally reacting (i.e., acoustic waves 
propagate through it normal to the surface), and the normalized acoustic impedance of the wall lining (generally a 
function of wall location) will be denoted by L. Throughout this study, all impedances are normalized with respect 
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Fig. 2 Schematic of inlet and aft-duct propagation and termination boundary, T f . 


It) the characteristic impedance of the air in the duel. The problem al hand is to predict the sound attenuation in the 
presence of the wall lining and grazing flow, given the mean flow field, sound source, normalized wall impedance, mid 
the duct termination condition. 

Codes available for predicting the attenuation in a 3D duct, such as is illustrated in Fig. 1, have generally been 
divided into non- vortical and vortical duct propagation codes. In this study, a suite of seven codes (four non- vortical 
and three vortical) is selected. The primary basis for code selection was that a two-dimensional (2D) and 3D version 
of the code be available. This decision was made so that the increase in wall clock time in transitioning from a 2D 
version to a 3D version of the same code could be accurately recorded. Although the Ac trail code 2 has a 2D and 
3D version, it is a commercial code and our current license agreement forbids our version to be exported to NASA’s 
Columbia Supercomputer. For this reason, Actran has not been included in the suite of selected codes. 


Ill, Non-vortical Duct Propagation Codes 

T Hh following section gives a brief overview of the basic assumptions, governing differential equations, boundary 
conditions, and numerical methods underlying development of the non- vortical duct propagation codes. 


A. Assumptions, Governing Equation, and Boundary Conditions 

The basic assumptions underlying the non- vortical duct propagation codes are (1) the acoustic field in the duct is 
linear, (2) the fluid in the duct is invisdd and irrotational, (3) the fluid is an ideal gas, and (4) the acoustic pressure and 
density are homentropicaUy related. These assumptions lead to a scalar wave equation (i.e., a converted Helmholtz 
equation) on the acoustic velocity potential, 1 2 7 


D! 




( 2 ) 


w here a time dependenee of the form e IW ? has been assumed, Cq is the sound speed, po is the mean static density, • 
denotes a dot product, is the irrotational mean velocity vector, is the gradient operator, / is the unit imaginary 
number, mid the material derivative is expressed in the frequency domain as 


^ = ko + 4^0 • ^ 


( 3 ) 


Upon obtaining the acoustic velocity potential, it is post-processed to obtain the acoustic pressure, /?, and acoustic 
particle velocity field, it 


P = 


Dh 

- po DC 


it = V<|) 


( 4 ) 
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The scalar wave equation ]Eq. (2)] requires three types of boundary conditions. One boundary condition is required at 
the source boundary, a second is required at the duel termination, and a third is needed along all rigid or acoustically 
lined wall segments. 

The source plane boundary condition is 

P = Ps, on r s (5) 

where p s is the known source plane acoustic pressure. At the duct termination, the boundary condition is specified in 
the form 

{p} = Kf] {pQColt •it}, along Ft (6) 

where [y is the dimensionless node impedance matrix 8 at points along T,, it is the outward unit normal vector to 
Ff, and {p}, { it } are vectors containing the acoustic pressure and the normal component of acoustic particle 
velocity at nodes along T/. Elements in the node impedance matrix, [y, may be chosen so that the termination is 
reflecting or nonreflecting. Further, the choice for elements in [y will be discussed in more detail in the results 
section of this report. 

There has been considerable discussion in the literature during the past half century concerning the acoustic bound- 
ary condition at an acoustically treated duct wall with a general shape mid mean flow'. However, the Myers boundary 
condition 9 has been generally accepted as correct in modem duct acoustic literature. Therefore, we choose to only 
include codes that implement the w all impedance boundary condition developed by Myers ? When written in terms of 
the primitive variables (i.e., the acoustic pressure and the acoustic particle velocity vector), this boundary condition is 

-p 0 O)TJW = f^-r)+l f 7t*7t*^lt 0 = Q on T w (7) 

Equation (7) is restricted to locally -reacting acoustic liners, and will not be valid for those liners that allow' wave 
propagation through the material in directions parallel to the flow' path. It is also valid along a rigid wall segment of 
duct, provided the normalized wall impedance is set to infinity (i.e., £ = oo). Note that the sound source [Eq. (5)], duct 
termination [Eq, (6)], and wall impedance [Eq. (7)] boundary conditions have been written in terms of the primitive 
variables so that they are consistent with the vortical flow’ boundary conditions that are presented in a later section. 
However, these boundary conditions [Eqs. (5) -(7)] are easily transformed to boundary conditions on the acoustic 
velocity potential by substitution from Eq. (4). 

B ♦ N urneri cal Methods 

It remains a formidable task to obtain solutions to the converted Helmholtz equation [Eq, (2)], coupled w ith the sound 
source, duct termi nation, mid wall impedance boundary conditions. Exact solutions to this equation for arbitrary 
boundary condition parameters, mid means flows have thus far not been found. However, the invention of large-scale 
digital computers has allowed numerical solutions to be obtained. These numerical solutions arc essential for optimum 
design of noise efficient aircraft. Over the past two decades, two methods for solving the 3D scalar wave equation (the 
converted Helmholtz equation), coupled to the sound source, duel termination, mid wall impedance boundary condi- 
tions, have emerged, 'iliese two methods are the parabolic approximation and the finite element method. Dougherty 1 
originally proposed the parabolic approximation for duct aeroacoustics and Watson 1 1 initially proposed the finite ele- 
ment methodology presented in this study. Although the original application of the finite element methodology 1 1 was 
for 2D impedance eduction in uniform flows, the method has recently been extended to 3D rectangular geometries and 
non-uniform mean flows, 5 The purpose of the follow ing subsections is to discuss these two numerical methods (i.e., 
the parabolic approximation and the finite element method) along with their salient features mid limitations. Read- 
ers are reminded that the details presented in this mid subsequent sections are provided as a quick review mid more 
technical detail is available in the original publications. lp 5f 6 

L The Paraho l ic Approxi matio n 

Recently, a Ian noise duct propagation and radiation code 12 was developed that utilizes a parabolic approximation 
to the scalar wave equation defined by Eq. (2). The code is composed of five distinct modules: (1) input/output 
specification, (2) fluid dynamics /acoustic grid generation, (3) background flow calculation, (4) acoustic propagation, 
mid (5) acoustic radiation. The acoustic propagation module is the focus of this section, mid it utilizes a parabolic 
approximation to the scalar wave equation. 1, 13 This approach aff ords very ef ficient propagation calculations, allowing 
for complex 3D geometries to be handled with relatively low computational costs. This efficiency comes at the 


4 of 20 



expense of reduced accuracy as the direction of propagation of an acoustic mode diverges from the preferred angle of 
the parabolic approximation. Additionally, loss of accuracy may occur when reflection of acoustic waves in the axial 
direction becomes important, as these effects are not captured in this formulation. Nevertheless, if appropriate care is 
taken to account for these limitations, the parabolic approximation can provide an ef ficient mid useful environment to 
perform 3D computations. 

Although the parabolic approximation mid general splitting technique were discussed in detail in previous work, 5 
a brief review is provided for the sake of completeness. The full formulation is based on a solution to Eq. (2) in 3D 
geometries. However, for simplicity, a rectangular duct is considered in which a uniform axial mean flow is present. 
In this situation, Eq. (2) is rewritten in the form 


-( 1 - JWg) H = V 2 J> - 2ikMo ^ + k 2 < |> 
dz A oz 


( 8 ) 


where V±,Mo, and k are the transverse Laplacian, mean flow Mach number, and freespace wave number, respectively 


d 2 d 2 Uq , 

Vl — Vj_ * Vj^ = — T — , Mo——, k— — 

dx 1 dy z Q} 


( 9 ) 


Here, wq is the uniform axial flow velocity mid an e~ iwt lime dependence is assumed. As discussed by Dougherty 13 
mid outlined by Corones, 14 a derivation based on splitting the solution of Eq. (2) into ‘positive’ traveling (i.e., <j) + ) mid 
‘negative’ traveling (i.e., (j)”) waves is proposed. This general splitting is represented as 


<J> + 

<T 


= [T] 



dZ 


= [T] 




♦ + 

<r 


( 10 ) 


where the splitting matrix, [T], is arbitrary. Substitution of the general variable splitting given by Eq. (10) into Eq. (8) 
leads to a coupled pair of differential equations for of the form 



Dougherty 1 considered plane waves to produce a splitting of the form 


r * 1 
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1 
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— ik 
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ik 

l-Mo J 

i 

l * J 


(ID 


(12) 


Typically, at this point in the derivation, the terms that couple the equations for <() ± are neglected to produce a single 
parabolic equation for the ‘positive 1 or ‘negative’ traveling wave. 

The general formulation outlined by Dougherty considers a system of orthogonal, curvilinear coordinates (^i , ^2, ^3) , 
defined so that ^3 is aligned with the steady potential mean flow. Following the aforementioned splitting process and 
neglecting the off-diagonal lenns in the resulting matrix equation gives the parabolic approximation 


alT 


I 13 ik 
1 +M 0 


- — vi^ 
2 ik x 


1 d 

2*, *2 Sfy , 


(hih 2 )$ + 


(13) 


where A, is the element of length given by 




2 

for i — 1,2,3. 


The source boundary condition is given by 

<|> + = <|>J, on r s (14) 

where, (j) ? is the known source potential. The solution for t|) 4 is marched from the source to the duel termination mid 
no termination boundary condition is required, as the ^ coupling is neglected. 


5 of 20 



The wall impedance boundary condition is also only approximately satisfied. This is achieved by assuming that 
the acoustic pressure and the normal component of acoustic particle velocity are continuous across an infinitely thin 
boundary layer. Denoting quantities within the duet interior with subscript I mid those below the infinitely thin 
boundary layer with subscript 2, the impedance boundary condition below the boundary layer may be written as 

+_L^! _ <$£ 

hj l 

where ^ * s the acoustic velocity potential below the boundary layer and the sign is taken as negative for minimum %j 
and positive for maximum %j. The continuity of pressure condition leads to 

Z)d>+ , D r i — y. 

= — ittxj>2 i =- = [-ico + lfo* v] 

where T^o is the mean flow velocity and the normal particle displacement condition gives 

(§) 


Finally, eliminating and neglecting derivatives of the Mach number in the direction of the flow leads to 


dty* ik hj , hjMo d /l\ , 

j K{\+m 0 ) 2 i+JMbafe Vt/ T 


(15) 


The parabolic equation |Fq. (13)] and associated boundary conditions [Hqs. (14) and (15)] can be efficiently solved 
using a marching method that advances the solution in the ^ direction. The classical fourth order Runge-Kutta 
method is well suited to this purpose and is used because of its accuracy and ease of implementation. T he transverse 
derivatives can be evaluated by any number of numerical schemes (e.g., finite dif ference, finite element). In the current 
implementation, a pseudospectral approach is used, in which the values of <|) + are approximated at the grid points by 
a series of sine mid cosine functions. 

IjcL mid ^ 2 be discretized on uniformly- Spaced grids taken from 0 to it. To approximate the transverse derivative 
of <|> + with respect to , the wave is modeled at each discrete value by 


tfi-1 


<fr + (li>^)=A 1 (| 2 )sm(l 1 )+A 2 {^)sin(2i 1 )+ J ft, (§ 2 ) cos (/i&) 

h=o 


(16) 


w r here Ai, A 2 , and #/,are modal coefficients and AT is the maximum grid index number in the direction. The 
w r all impedance boundary condition from Fq. (15) is used to determine the coefficients Aiand A 2> and a fast Fourier 
transform method is used to evaluate the coefficients. Once these coef ficients are determined, the derivatives are 
calculated from 

= A, (1 2 )cos(|i) + 2A 2 (| 2 )cos(21 1 )- 2 hB h (^)sin (Ifo) (17) 

6= 0 


mid 


-£r = ~ A ' (l2)sin(ii)-4 J 4 2 (| 2 )sin(2|i)- £ l\B h (^cos^lO 


/,-0 


(18) 


T he same process is followed to evaluate the derivatives in the % 2 direction. Substituting the transverse derivatives into 
Fq, (13) completes the evaluation of the right-hand side and the solution is marched forward to obtain <j) + . 


2. The Finite Element Method 


This section describes the weak, finite element solution to the convected Helmholtz equation [Fq. (2) [ that satisfies 
the source, duct termination, and wall impedance boundary condition. To begin, the weak, Galerkin, finite element 
approximate solution to the acoustic potential, <j>, requires that the acoustic potential he first approximated by a class 
of complete basis functions, Nf, 

NN 



( 19 ) 


6 of 20 



where <E»/ is a set of generalized coordinates whose relationship to § will be determined and the basis functions are 
understood to be function of the spatial coordinates, (z, Jt,y). The residual error, R\ , resulting from the substitution of 
Eq. (19) into the governing differential equation [Fiq. (2)] is now minimized, by requiring that it be orthogonal to each 
of the basis functions, Nj 


and the residual error is 


L 


RiNjdQ = 0 


Ri = - 


D 

Dt 



( 20 ) 


( 21 ) 


Equation (20) is a system of NN linear equations in NN unknown generalized coordinates, d>/. The weak formulation 
is achieved by applying Green’s theorem to the second derivative terms in Eq, (20), and reducing them to first order. 
Eor example. 


jW • (poV+) NjdQ =~ J Po (V«|* • Tf ) (^Nj • it'jdQ + J p 0 % • T?NjdT 
Recalling that, it • it = • it [see Eq. (4)], gives 


( 22 ) 


J | 7 * (po^) NjdQ = - J po (^<t> • (^Nj • ^ dQ + j p^lt • itNjdr 


(23) 


Substituting the termination [Fq. (6)] and wall impedance [Fq. (7)] boundary conditions into the surface integral in 
Fq. (23), followed by substitution of Eq. (23) into Fq. (20), allows Eq. (20) to be written in weak form. 

The system of AW linear equations in AW unknowns resulting from the weak form of Eq. (20) is constructed using 
a 3D finite element methodology. The finite element methodology used in this paper guarantees continuity of the 
acoustic potential and the acoustic particle velocity vector (i.e., V(|)) within and on the boundaries of £2 The 3D finite 
element method assumes that there are N grid lines in the axial or z direction, M grid lines in the x direction, mid Q 
grid lines in the y direction of the duct. As shown in Fig 3., grid lines need not be equally spaced mid an element is 
identified by a triple index notation, [/,/, £]. A typical element, [7,7,W], is bounded on the right and left by the grid 
lines, y = y^, and, y = y ^+ 1 , on the bottom and top by grid lines, x = xj, and, x = xj+\ , and on the front and back by 
grid lines, z = Zu and, z — Z/+ 1 - Brick elements with eight nodes are used in the analysis. A typical brick element is 



Fig. 3 Finite- element discretization of 3D duct. 

shown in Fig 4. Each brick element, [7,7, Af], consists of eight local node numbers labeled 1, 2, ... 8. Each element 
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Fig* 4 A Typical 3D finite- el erne lit brick element and local node numbering system* 


is considered to have length a , height b , and width d , as shown. The objective of the 3D finite element method is to 
obtain the unknown acoustic velocity potential and its derivatives at the nodes of each of the (N — 1 ) (M — 1 ) (Q — 1 ) 
elements. Within each element, 4> is represented as linear combinations of sixty-four, 3D, cubic, Hermite polynomials, 
N \ , N 2 >N$ - ** ^64 that comprise a complete set of basis functions within each element 


1=64 

<t»= ^ (24) 

The finite element analysis allows the sixty-four 3D basis functions, Af/, to be easily constructed from the products of 
the four 1 D cubic Hennite polynomial basis f unctions 


= 1 -3(z/a) 2 + 2(z/a) 3 , Ji{z,a) = z{z/a- l) 2 
Mz,a) = {z/a) 2 (3 -2zfa), f*{z,a) = (z 2 /a){z/a- 1) 

in the three coordinate directions. For example, it is easily shown that 

N { =A(z,a)Mx,b)My,d), N 2 = f 2 (z,a)M \x,b)My,d), ... N 64 = f 4 {z,a)f 4 (x,b)f 4 (y,d) (26) 


Similarly, the relationship between the unknown generalized coordinates, d>/, and the velocity potential at each node 
points, [zuxj^yK), is easily obtained. For example, 


<J>1 = Jk), 


d?2 — 


W(zi,xj,yic) 

dz 


*&64 = 


dzdxdy 


(27) 


Substituting Fqs. (24) into Eq. (20), using Green’s theorem to reduce the second order derivatives to first order 
derivatives, and then substituting the wall impedance and duet tenni nation boundary conditions into the surface in- 
tegrals, allows one to compute the contributions to the minimization of the total residual error from each individual 
element 


W-lAZ-lO-l 

j RiNjdQ = V V V [A I J - K ]{q>\!dA} (28) 

where the sum on the right-hand side of Hq. (28) is performed using the usual rules of finite element assembly. In 
Fq (28), [aIM.^ 1] is a 64x 64 complex matrix for each element, [/,/, and } is a 64 x 1 column vector 

containing the unknown acoustic potential and its derivatives at the eight nodes of the element. The coef ficients in the 
local stiffness matrix, are computed in closed form using the Mathemalica intrinsic functions. 15 
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Assembly of the global equations for the computational domain is a basic procedure in the finite- element method. 
Appropriate shil ling of rows mid columns is all that is required to add the local element matrix, [A^’^], directly into 
ihe global matrix, [Af]* 16 This process leads to a set of linear matrix equations of the form 

MW = {0} (29) 

where [Ay] is a complex, block tri diagonal matrix whose order is SMNQ, mid {d>} and {0} are 8 MNQx 1 column 
vectors. The vector (4>} contains the nodal values of the unknown acoustic potential mid its derivatives at the nodes 
of the elements and {0} is the null vector, llie matrix [A/] is singular until the source conditions are applied. 

The source is introduced into Eq, (29) by constraining the nodal degrees of freedom at the source plane in the usual 
manner* 16 After implementation of the source condition, Eq. (29) remains block tridiagonal and is of the form 

[A]{$} = {F} (30) 

where {F} now contains the source effects and [A] is a modified form of [A/] with the same structure. Much practical 
importance arises from the matrix structure in Eq. (30), as it is convenient for minimizing storage and maximizing 
computational efficiency* Special matrix techniques exist for a solution of this structure. An imsymmetri cal parallel 
sp;irse solver with equation reordering to minimize fill is used to decompose [A] into the products of a unit lower 
triangular, a unit upper triangular matrix, and a nonsingular diagonal matrix. The sequential operations of forward 
and backward substitution are used to obtain the solution vector, { } . All computations are performed on a parallel 
computer and only the nonzero elements in [A] are stored and operated upon during the solution stage. 

IV, Vortical Duct Propagation Codes 

T HE following section gives a brief overview of the basic assumptions, governing diffcrcnlkil equations, boundary 
conditions, and numerical methods underlying development of the vortical duet propagation codes. 

A. Assumptions, Governing Equations, and Boundary Conditions 

The vortical flow duct propagation codes are based on the same assumption as that for non -vortical flow, except that 
the irrotationality assumption is not imposed. Consequently, there is not mi acoustic velocity potential mid the solution 
is Obtained in terms of the primitive variables that satisfy the sound source, wall impedance, mid duet termination 
boundary conditions* The basic assumptions underlying code development are that (1) the acoustic field is linear, 
(2) the fluid flowing in the duct is inviscid, (3) the fluid is an ideal gas, and (4) the acoustic pressure and density 
remain homentropically related. 

Using the above assumptions, the linearized equation that governs conservation of mass in a flow duct with ail 
inviscid fluid is 3,4 

(?) + + = ° ( 31 ) 

and the linearized equations that govern the conservation of momentum i s 

ypo^F-+ypo(jt»^)^o+p^^-+(o^p = o (32) 

where p 0 is the mean static pressure, and y is the ratio of the specific heat at constant pressure to that at constant 
volume, 'Hie acoustic particle velocity vector, it, is assumed to contain components u , v, and w in the x 7 y, and z 
directions, respectively. In 3D, the linearized continuity equation [Eqs. (31)] is a scalar equation and Eq. (32) is a 
vector equation with three components. Thus, the linearized mass mid momentum equations are a system or four 
equations in four unknowns (i.e., p, h ? v, and w) that can be solved uniquely to determine the unknow n variables that 
satisfy the sound source, duet termination, and wall impedance boundary conditions. 

In the current subsonic flow problem, the linearized aeouslie system in three space dimensions requires three 
boundary conditions ai the source boundary, one boundary condition at the duet termination, mid another boundary 
condition at either a rigid or acoustically lined wall segment. 17 The wall impedance boundary condition remains the 
Myers boundary condition [Eq. (7)], and we w ill specify the acoustic pressure and the transverse components of the 
acoustic particle velocity at the source boundary as 
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where p s , u s and v s are the known values of p, u and v at points along r s . Finally, at the duct termination, the boundary 
condition remains that given by Fq. ( 6 ). 


B . N umeri cal Methods 

Four vortical flow, duct propagation codes, are included in the current suite. Bach code solves Fqs, (31) and (32), and 
satisfies the sound source [Eqs. (33)], duct termination [Fq. ( 6 )], and wall impedance [Fqs. (7)] boundary conditions. 
The first two codes arc based upon a linite element solution methodology, and were initially developed at the NASA 
Langley Research Center in 2006. 4 The second two codes use a pseudo -time marching approach, mid were developed 
at the Boeing aircraf t company with partial funding from NASA in 2006. 3 The following subsections give a brief’ 
description of the numerical methods underlying the development of the vortical flow', duel propagation codes. 


1. The Finite Element Method 

The two vortical, finite element, duct propagation codes solve Fqs. (31) and (32) using a linite element methodology. 
The numerical procedure parallels that used in section 2. However, only linear basis functions are used mid the wall 
mid exit impedance boundary conditions are satisfied by constraining the nodal degrees of freedom. To begin, the 
Galerkin, linite clement approximate solution to the linearized mass mid momentum equations, requires that each 
primitive variable be first approximated by a class of complete basis functions, Nj, 

NN NN AW AW 

p = ^N/pi, u = ^N[U[, v = ^ Njv/, w = ^ N/Wj (34) 

where and >v/ are generalized coordinates. The residual errors, R r mid {/? WI } resulting from the substitution of 

Fq. (34) into the linearized mass and momentum equations [Fqs. (31) mid (32)] arc now ? minimized, by requiring that 
they be orthogonal to each of the basis functions, Nj 

[ RcNjdQ = 0 (35) 

J&2 


f {R m }NjdQ = Q 

Jii 

Rc = CqT (p'j + YPo ^ •Tt + p$»Tt 0 + 4y(t^*^) ^ 


(36) 

(37) 


[ R m} = jypo^ + YPo ^o+p^^+co^pj 

The rcsidiud error, R Ci is a scalar, while the residual error, {/?^} , is a vector (with three components) 


( 38 ) 


R* 

Rn 


(39) 


where Rmi^Rmx ^d Rmy arc the residual error components in the z, x mid y directions, respectively, of the linearized mo- 
mentum equation. Equations (35) and (36) reduce to a system of 4AW linear equations in 4 AW unknown generalized 
coordinates, pj , uj , v/ and wj . 

The system of 4 AW linear equations mid 4AW unknown generalized coordinates is constructed using a 3D linite 
element method similar to that used for the convected Helmholtz equation. The 3D finite element method assumes that 
there are N nodes in the axial or ; direction, M nodes in the x direction, and Q nodes in the y direction of the duct, just 
as in Fig. 3, Brick elements with eight nodes are again used (Fig. 4) and the objective of the 3D finite element method 
is to obtain the solution for the four unknown primitive variables at the nodes of each of the (N — l)(M — 1 ) (C? — 1 ) 
elements. Within each element, each primitive variable is represented as a linear combinations of eight, 3D, linear 
polynomials, N\ W 2 W 3 . . Wg that comprise a complete set of basis functions within each element. Just as with the 
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scalar wave equation, the finite element method allows the eight 3D basis functions, Afy, and generalized coordinates 
Ph u U v /, and w/ to be easily constructed from products of the two linear polynomial basis functions 


h (Z, a ) = 1 “ z/a, f 2 (z, a)=z/a (40) 

in each of the three coordinate directions, along with the requirement that the primitive variables each be continuous 
within Q and along inter-element boundaries. For example, it can be shown that 

Ni = fx (z, a)/, (*, b)h (y, d ) , N 2 = f\ ( Zt a)fi (x y b)f 2 (y y d) t ... N s = f 2 (z y a) ft (*, b)f x (y, d) (41) 

and that each generalized coordinate, and wj are therefore the values of p,u 7 v and w at local node L Sub- 

stituting the series expansions in Fq. (34) into Eqs. (35) and (36) allows one to compute the contributions to the 
minimization of the residual from each individual dement 


/ 


Rc 

Rmz 

Rmx 




) NjdQ = 


J 


N— l M— t Q— 1 


(42) 


where the sum on the right-hand side of Eq. (42) is again performed using the usual rules of finite element assembly. 
In Eq. (42), [A^ 1 ^] is a 32x32 complex matrix for each element, [7,7,^], and {cpIVAl} is a 32 x 1 column vector 
containing the unknown primitive variables at the eight nodes of the element. The coefficients in the local stiffness 
matrix, are again computed in closed form using the Mathematica intrinsic functions, 15 Assembling the 

elements for the entire domain results in a matrix equation of the form: 


= { 0 } 


(43) 


where [Ay] is a complex, block tri diagonal matrix, whose order is 4MNQ , and {d>} and {0} arc 4MNQ x 1 column 
vectors. The vector {4>} contains the nodal values of the primitive variables at the nodes of the elements and {0} is 
the null vector. 

Before a non-trivial solution to Eq. (43) can be obtained, source, termination plane, and wall impedance boundary 
conditions must be imposed. To satisfy the source, wall impedance, and termination plane boundary conditions, the 
residual error corresponding to these boundary conditions is minimized [i,e., after substitution from Eq. (34)], For 
example, to satisfy the source boundary condition, the source residual, R St is minimized 



(44) 


R s = 


P~Ps 

u-u s 


(45) 


{ v-v s ) 

Equations (44) constitute a set of constraints that are imposed on the solution to Eq. (43). Similarly, to satisfy the wall 
impedance boundary condition, we minimize the wall residual, R w 



(46) 


Equations (46) constitute a set of constraints that, when imposed on the solution to Eq. (43), satisfy the wall impedance 
boundary condition. A similar procedure is used to satisfy the te rmina tion boundary condition. The above methodol- 
ogy of introducing the boundary conditions as constraints results in a modified set of matrix equations of the form 


[A]m = {F} (48) 

where {F} contains source cffeels and [A] has a structure similar to [Ay]. Equation (48) is also solved using a parallel 
band solver to obtain the solution vector, {<b}. 
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2 . The Pseudo-Time Method 


A pseudo-time method has also been developed for solutions of Eqs. (31) and (32) when coupled with the source 
[Eq. (33)], the duet termination [Eq. (6)], mid wall impedance [Eq, (7)] boundary conditions. The starting point for 
the pseudo -time marching code or method is to rewrite Eq. (31) mid Eq. (32) in a compact vector format at a location 


-^{Q}+[A]{~ } + [s]{^} + [q{f }+PK0 = o, \Q) = 


m 

u 

V 

w ) 


(49) 


where [A] > [5] and [C] are 4 x 4 matrices that contain mean flow effects and [£)] is a 4 x 4 matrix containing terms 
related to the gradient of the mean flow field, The coefficients in [A], [B] 7 [C], and [D] are rather lengthy and are not 
written explicitly in this report. Detailed expressions for these matrices can be found in the original reference. 3 

Adding a pseudo-time derivative to the compact vector form of the linearized equations that govern continuity of 
mass and momentum |Eq. (49)] gives 


{ = -»»{£} + [AH ^ } + [B]{ + [C]{ + [D]{Q} (50) 

where t is a pseudo-time variable. The pseudo-time derivative is introduced so that Eq, (50) may be marched in pseudo- 
time until the solution is steady. Thus, the steady solution to Eqs. (50) is the solution to Eq. (49), Equation (50) is 
discretized in pseudo -time using a forward- Euler scheme 


r C n+ 1 -g n 
1 At 


} + /co{G M+1 } = [A]{ + [B]{ + [D]{QT} 

ox dy dz 


( 51 ) 


where n is the pseudo- lime index. At is the pseudo-time step, and the lirst term on the right side of Eq. (50) is treated 
implicitly. Equation (51) is now written in a delta form 


{AG} 


At 


+ to] = -to{G”> +[ a]{^}+ + [c]{^} + mm, a q = g w+i 


-or 


(52) 


Upon approximating the spatial derivatives and then implementing the boundary conditions, Eq, (52) eau be marched in 
pseudo-time until the solution is steady (i.c., until {A(9} = {0}) to achieve the frequency domain solution. The spatial 
gradients are approximated using a fourth order optimized finite difference technique, referred to as the dispersion- 
re! alien -preserving (DRP) scheme. 18 The DRP scheme is a central difference approximation to the gradient in each 
coordinate direction that uses a seven-point stencil. The dispersion relation- preserving scheme is especially designed 
to minimize dissipation and dispersion errors with a resolution of eight points per wavelength. At and near boundaries 
of the computational domain, the DRP central difference scheme requires points from outside the computational 
domain, This is avoided by switching to a fourth order, optimized, one-sided difference at and near these boundaries. 
Further, an eighth-order filtering scheme is applied to prevent spurious waves from contaminating the solution domain 
as the solutionis advanced forward m pseudo- time, The sound source, wall impedance, and duct termination boundary 
conditions are inserted into the pseudo-time marching method in strong form. Other intricate details involved in the 
implementation of the DRP scheme to the linearized equations of mass and momentum continuity, as well as the 
implementation of other acoustic boundary conditions not considered in this paper arc described in detail in another 
reference. 3 

The pseudo -time marching method is implemented on a curvilinear mesh using uniform computational coordi- 
nates. Thus, solutions using arbitrary geometries are possible. In addition, the pseudo -time marching code uses a 
domain decomposition scheme that utilizes the message passing interface (MPI) method to reduce wall clock time on 
parallel computers. The domain decomposition is only in the primary flow direction; thus, the other two dimensions 
are not parallelized. Therefore, the 2D version of the pseudo-time marching method is not currently parallelized. 


V. Liner Attenuation 

U PON obtaining the solution vector, {4>}, it is post-processed to obtain the acoustic velocity potential in non- 
vortieal flows and the primitive variables (i.e,, p, u, v, and w) in vortical flow s. Because only the acoustic pressure 
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is readily measurable, it is convenient to use the measured difference in the integrated sound pressure level between 
the source plane tmd an axial plane , z = Zi, to compute the liner attenuation 


mpl=iou, hJi£) 

(53) 

SPL s = f pp*dT 
J r. 

(54) 

SPL zl = [ pp*dT 

(55) 


where the superscript, *, denote the complex conjugate. Equation (53) is used as a measure of liner attenuation 
throughout this report. The location of the axial plane, z — Zl [Eq. (55)], is taken either upstream of the source plane 
for inlet propagation or downstream of the source plane for aft propagation. 

VI* Results and Discussion 

S EVEN duel propagation codes (3 non- vortical and 4 vortical) were exported to the Columbia clusters to perform 
ihe current assessment. Table 1 gives an overview of the code names (CH2DDS, LEE2DDS, LEE2DIS, CH3DPA, 


Table 1 . Summary of duct propagation codes exported to the Columbia clusters. 


Name 

Dimensionality 

Governing Differential Equation 

Numerical Method 

CH2DDS 

2D 

Convected Helmholtz 

Cubic Finite Element 

LEE2DDS 

2D 

linearized Euler 

Linear Finite Element 

LEE2DIS 

2D 

linearized Euler 

Pseudo Time- Stepping 

CH3DPA 

3D 

Convected Helmholtz 

Parabolic Approximation 

CH3DDS 

3D 

Convected Helmholtz 

Cubic Finite Element 

LEE3DDS 

3D 

linearized Euler 

linear Finite Element 

LEE3DIS 

3D 

linearized Euler 

Pseudo Time-Stepping 


CH3DDS, LEE3DDS, LEE3DIS), as well as the dimensionality (2D or 3D), governing differential equation (con- 
veeled Helmholtz or linearized Euler), and numerical method (parabolic approximation, finite element, or pseudo- 
time domain) underlying the development of each code. Note that the parabolic approximate solution to the convected 
Helmholtz equation is run only in 3D to limit the number of computer runs and because the 3D and 2D code run at 
nearly equal speeds. Further, the duct propagation code consisting of a 3D finite element solution to the linearized 
Euler was not run during this assessment due to time constraints. 

Each code is run for a geometry corresponding to the LaRC Crazing Incidence lube (GIT) 11 using an assumed 
uniform mean flow. This geometry and mean flow was chosen because: (1) approximate liner attenuations from 
mode theory are available for comparisons, (2) experimental data are available to compare with code predictions, mid 
(3) some of ihe duel propagation codes are currently written for a rectangular geometry such as present in the GIT. The 
GIT has a square cross-section, with a height and width of 0.051 m. For this study, the distance between source mid 
measurement planes is 0.61 m (i.e., <^=0.61 m). Sound mid flow in the GIT are currently directed in the same direction, 
si) that only aft -duet propagation was modeled in the duet propagation codes. A grid refinement study showed that 
Jill codes were nearly converged on a uniform computational grid with 97 points in the axial and 21 points in the two 
transverse directions for the Mach numbers and frequencies used in this study. Thus, all 3D and 2D codes were run on 
uniform 97 x 2 1 x 2 1 and 97 x 21 baseline grids, respectively. 

A, Accuracy Assessment 

The most important metric in assessing a propagation code is the accuracy of the predicted attenuation. To assess 
the accuracy of the predicted attenuation, two example problems were chosen for which the attenuation could be 
determined analytically. The first example problem tests the ability of the duel propagation codes to predict the 
attenuation of a plane wave propagating through a hardwall duct with flow'. It is w ell known that a plane wave does not 
attenuate as it propagates along a hard wall duel, mid therefore the codes should predict zero attenuation. The second 


13 of 20 




example problem tests the ability of each code to predict the finite attenuation of a cut -on mode in a soft wall duct with 
a uniform lining. For this case, an approximate value for the attenuation is obtained using mode theory 

Table 2. Attenuation in dB of plane wave in a hardwall duct with an aneclioic termination 


Mo = 0.0 


/.kHz 

Exact 

C1I2DDS 

LEE2DDS 

LEE2DIS 

CH3DPA 

CII3DDS 

LEE3DIS 

1.0 

0.00 

0.00 

0.01 

0.00 

0.00 

0,00 

0,00 

2.0 

0.00 

0.00 

0.02 

0.00 

0.00 

0.00 

0.00 

3.0 

0.00 

0.00 

0.00 

0.05 

0.00 

0.00 

0.09 


M 0 = 03 


1.0 

0.00 

0.00 

0.01 

0.00 

0.00 

0,00 

0.00 

2.0 

0.00 

0,00 

0,01 

0.00 

0.00 

0,00 

0,00 

3.0 

0.00 

0,00 

0,21 

0.01 

0.00 

0,00 

0,01 


M 0 = 0.5 


1.0 

0.00 

0,00 

0,00 

0.00 

0.00 

0,00 

0,00 

2.0 

0.00 

0,00 

0,02 

0.00 

0.00 

0,00 

0,00 

3.0 

0.00 

0.00 

0.06 

0.00 

0.00 

0,00 

0.00 


In the first example problem, the hardwall attenuations arc simulated for nine diff erent conditions obtained f rom 
three selected excitation frequencies (1.0, 2.0, mid 3.0 kHz) mid three selected mean flow' Mach numbers (0.0, 03, 
0.5). For each simulation, the source is planar (i.e., p s = 1), the termination is lioiircflecliiig (i.e. [y = [/]), mid 
the niemi static pressure, density, and temperature are those at standard atmospheric conditions (po— 10 1325.0 N/m 2 , 
Po=l 3 kg/m 3 , Jb=293.0 K). Here, [/] is the identity matrix. Results are presented in Table 2, Computer simulations 
from each of the selected codes are in good agreement with the exact value of zero. There is a small attenuation of 
0.21 for LEE2DDS at Mach 0.3 for the highest frequency (i.e., 3.0 kHz). r fhe discrepancy in this case is apparently 
due to the use of a low -order approximation (i.e., expansion in terms of linear basis functions) to the acoustic pressure 
for the highest excitation frequency. When the LEE2DDS grid density w as increased by a factor of two, the predicted 
attenuation was reduced from 0.21 to 0.05 dB. 

The second example tests the ability of each code to predict the finite attenuation of a cut- on mode in a soft wall 
duet with a uniform lining [i.e., a duet for which the lower and two sidewalls are hard mid the upper wall contains 
a uniform liner (Impedance values for the softwall are chosen from a database of’ six production liners that were 
tested in the GTT. Impedance values were judiciously chosen from the database so that they covered the maximum mid 
minimum ranges of resistance and reactance available in the database of six liners)]. In this ease mi approximate value 
for the attenuation can be obtained using mode theory. It is easily shown that the exact value of the sound pressure 
level attenuation from Eq. (53) for this configuration is 

A SPL = 20iog lo [e] %[K nZl J (56) 

w'here 3[ ] denotes the imaginary part of the complex expression enclosed within the brackets, zl is the location 
of the downstream plane (i.e., 0.610 in for the current study), and the axial propagation constant K n satisfies the 
transcendental equation 

ika„ = tan (X„/f) (57) 

w here H is the duel height and 



The trans verse wavenumber, X„, and axial propagation constant, K n , are related by the dispersion relation 

K, «o±\/l — (1 -Mg)(V,/*) 2 

k - -(i-K) (39) 

The axial propagation constant, K ny is approximated by finding the roots of the transcendental equation [Fq. (57)] using 
a Newton-Raphson iterative scheme. Note that for a given value of X n , the axial propagation constant K n satisf ying 
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Eq, (59) has two roots. Here, we will choose the root that corresponds to right -moving waves in the duct. It is easily 
shown that right- moving waves are identified as those waves for which the axial propagation constant, K, u has a zero 
or negative imaginary part and a positive real part, respectively. In addition, the above solution [i.e. s Eq. (56)-(59)] 
considers only the zeroth order spanwise mode so that the solution in the third dimension (i.e., the y dimension) is 
planar. Thus, the solution for this softwall example problem is only 2D ( p — see Fig. 3), so that both the 2D 

and 3D version of a code should return identical solutions. 

Attenuations due to a softwall were computed for 36 different conditions, These 36 conditions correspond to 
identical frequencies and Mach numbers as in Table 2. However, four softwall impedances were used with each 
combination of Mach number and frequency, The sound source, p s , w as chosen to be the least attenuated mode in 
the duct, the termination w r as chosen to be nonreflecting (i.e. [y = E;^/]), and the mean static pressure, density, 
and temperature are those at the standard atmospheric conditions used in the hard wall simulations. Here, is the 
characteristic impedance for the least attenuated right- moving wave and was computed from the exact expression for 
the outgoing wave impedance 


Ke 


k-MpKn 

K* 


(60) 


Attenuations predicted from CH2DDS, LEE2DDS, LEE2DIS, CH3DPA and CH3DDS for the softwall duel with 
mi anechoic termination and zero flow are compared to the approximate attenuation from mode theory in Table 3* The 
normalized resistance and reactance of the wall lining are given in columns 2 and 3, respectively, of the table. Recoil 1 
that the sound source in this second example was a softwall duct mode. Results from LEE3DIS were not computed 
because the available version of LEE3DIS was designed to take source inputs in the form of hard wall duet modes 
only As shown in Table 3, CH2DDS, LEE2DDS, CH3DDS, and LEF2DLS arc in good agreement with mode theory 


Table 3* Attenuation in dB of lowest mode in a softwall duct with an anechoic termination 


M 0 = 0.0 


kHz. 

*K1 

m 

Theory 

CH2DDS 

LEE2DDS 

LEE2DIS 

CH3DPA 

CH3DDS 

1.0 

0.50 

-0.50 

61.70 

61.70 

59.50 

58.62 

68.77 

61.10 

1.0 

0.50 

0.50 

50.00 

50.00 

49.80 

48.70 

43.15 

49.50 

1.0 

2.00 

-0.50 

24.60 

24.60 

24.40 

24.89 

25.73 

24.30 

1.0 

2.00 

0.50 

24.00 

24.00 

23.90 

23.77 

24.25 

23.80 

2.0 

0.50 

-0.50 

143.00 

143.00 

61.70 

73.27 

79,00 

142.00 

2.0 

0.50 

0.50 

19.30 

19.30 

19.20 

18.13 

35.52 

19.10 

2.0 

2.00 

-0.50 

26.50 

26.50 

26.30 

26.60 

25.39 

26.20 

2.0 

2.00 

0.50 

21.00 

21.00 

20.90 

20.75 

22.71 

20.80 

3.0 

0.50 

-0.50 

36.10 

36.10 

35.90 

33.88 

24.06 

35.70 

3.0 

0.50 

0.50 

10.30 

10.30 

10.20 

9.76 

8.35 

10.20 

3.0 

2.00 

-0.50 

26.40 

26.40 

26.30 

25.87 

27,63 

26.10 

3.0 

2.00 

0.50 

17.90 

17.90 

17.70 

17.35 

17,63 

17.70 


results with one exception. Note that at 2.0 kHz, LEE2DDS and LEE2DIS fails to predict the large attenuation of 143 
dB when the impedance is £ = 0,5 — 0,5/. Further, the finite element codes (C112DDS, CII3DDS, and LEE2DDS) 
were observed to converge to the mode theory attenuation from below (i.e., the attenuation is slightly lower than 
that of mode theory on a coarse grid, but gradually increases to the mode theory result as the grid is refined). The 
failure of three of the codes (LEE2DDS, LEE2DIS mid CH3DPA) to predict the large peak attenuation of 143 dB was 
somewhat surprising and should be further investigated. LEE2DIS is observed to be slightly less accurate than the 
duel propagation codes that use a direct solve strategy (CH2DDS, CH3DDS, mid LEE2DDS), Given the low fidelity 
of CH3DPA compared to the other codes in the suite, its perf ormance at f requencies away from the 143 dB peak 
is good, except for the 2.0 kHz, results when ^ = 0,5 + G.5i mid the 3.0 kHz results when ^ = 0,5 — 0.5r, However, 
these two frequencies and impedances do correspond to the sources with the widest wave angles where the parabolic 
approximation is less accurate, 

Softw all attenuations predicted at Mach 0.3 and Mach 0.5 are given in Tables 4 and 5, respectively Several 
observations are noted for these Mach numbers. A significant reduction in the peak attenuation is observed at 2.0 kHz 


15 of 20 




Table 4 Attenuation in dB of lowest mode in a softwall duct with an anechoic termination* 


M 0 = 03 


kHz 

m 


Theory 

CH2DDS 

LEE2DDS 

LEE2DIS 

CH3DPA 

CH3DDS 

1.0 

0.50 

-0.50 

21.20 

21.20 

20.80 

26.90 

39.76 

21.00 

1.0 

0.50 

0.50 

44.10 

44.10 

44.10 

37.48 

29.85 

43.70 

1.0 

2.00 

-0.50 

13.70 

13.70 

13.60 

18.51 

15.57 

13.50 

1,0 

2.00 

0.50 

15.00 

15.00 

15.00 

18.79 

15.03 

14.90 

2.0 

0.50 

-0.50 

39.10 

39.10 

40.70 

66.93 

43,01 

38.70 

2.0 

0.50 

0.50 

19.20 

19.20 

19.30 

18.49 

24,90 

19.00 

2.0 

2.00 

-0.50 

15.20 

15.20 

15.10 

20.41 

15.31 

15.00 

2.0 

2.00 

0.50 

14.00 

14.00 

14.00 

17.25 

14.28 

13.80 

3.0 

0.50 

-0.50 

52.30 

52.30 

52.50 

43.14 

56.22 

51.70 

3.0 

0.50 

0.50 

11.50 

11.50 

11.70 

10.51 

9.78 

11.40 

3.0 

2.00 

-0.50 

15.80 

15.80 

15.80 

20.52 

16,73 

15.60 

3.0 

2.00 

0.50 

12.80 

12,80 

12.90 

15.04 

12,72 

12.70 


Table 5 Attenuation in dB of lowest mode in a softwall duct with an anechoic termination* 


M 0 = 0.5 


kHz 

m 


Theory 

CH2DDS 

LEE2DDS 

LEE2DIS 

CH3DPA 

CH3DDS 

1.0 

0.50 

-0.50 

13.30 

13.30 

13.00 

23.57 

31.12 

13.10 

1.0 

0.50 

0.50 

38.90 

39.00 

39.10 

41.33 

25.53 

38.50 

1.0 

2.00 

-0.50 

9.93 

9.93 

9.87 

15.78 

12,05 

9,83 

1.0 

2.00 

0.50 

11.40 

11.40 

11.40 

16.74 

11,74 

11.30 

2.0 

0.50 

-0.50 

23.20 

23.20 

23.10 

26.28 

31,47 

23.00 

2.0 

0.50 

0.50 

18.50 

18.50 

18.70 

18.85 

20,48 

18.30 

2.0 

2.00 

-0.50 

11.10 

11,10 

11.10 

17.70 

11.63 

11.00 

2.0 

2.00 

0.50 

10.90 

10.90 

10.90 

15.46 

11.03 

10.70 

3.0 

0.50 

-0.50 

46.80 

46.80 

45.60 

46.12 

51.67 

46.30 

3.0 

0.50 

0.50 

11.60 

11.60 

11.90 

10.81 

9.80 

11.40 

3.0 

2.00 

-0.50 

11.60 

11,60 

11.60 

17.81 

12,44 

11.50 

3.0 

2.00 

0.50 

10.20 

10.20 

10.30 

13.72 

10,12 

10.10 
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for t — 0,5 — 0.5; as the Mach number is increased (see lables 3-5). CH2DDS, LHH2DDS and CII3DDS remain in 
good agreement with the mode theory results (note that at 2.0 kHz that LEE2DDS is in good agreement with the mode 
theory results now that the peak attenuation is much lower). Also, given the low fidelity of CH3DPA eompzired to 
the other codes, its predicted attenuation is good, except at a few selected frequencies where wide-angle effects are 
more pronounced. Finally, the attenuation of LEIi2DIS is diverging from the attenuation predicted by mode theory 
and from the attenuation predicted by other codes at these Mach numbers. This is mi unexpected result but appears due 
to LEE2DIS nonreflecting boundary condition that is less accurate for the nonplanar wave source used in this second 
example problem. Note that LEE2DIS and LEE3D1S are in good agreement with exact results for the planar source 
example presented in 'fable 1. 

Our first example problem contained a hardwall and the second example contained a uniform liner. In these two 
example problems, the inputs to the prediction codes w ere known and exact or approximate solutions for the attenua- 
tions were determined analytically for comparison with code predictions. A third example is now presented for which 
the lower wall, two sidewalls, and portions of the upper wall are rigid, A 0.41 m-Iong, uniformly treated section is 
positioned 0. 10 m downstream of the source plane. The lest liner was a conventional perforate- over-honey comb ma- 
terial whose grazing flow 7 impedance was obtained from measurements in the GIT This configuration is referenced in 
the following discussion as the segmented liner example. To avoid reflections from the trailing edge of the segmented 
liner, the test section length was increased to 0.71 1 m and the grid spacing was kepi at that of the previous examples 
by increasing the grid density to a 1 13 x 21 evenly spaced grid. Unfortunately, approximate analytical expressions 
for the segmented liner attenuation are not av ail able. Predicted attenuations are therefore compared to measured data. 
Additionally, because the code inputs must be measured, they are not known exactly but are subject to measurement 
uncertainty. Uncertainty also exists in the measured attenuation. Therefore, measurement and predictions are com- 
pared on a statistical basis. 

For the segmented liner, attenuation predictions are performed for a 140 dR plane wave source in the GIT, in which 
the upper wall contains a 0.40 m-long liner. Results were acquired at three flow' rates for frequencies of 1.0, 2,0 mid 
3.0 kHz. Data were acquired at each of these conditions at least eight times in order to allow for a small -sample 
statistical evaluation of code inputs and measured attenuation. These tests w'erc conducted at different times over a 
period of one year so as to minimize systematic errors due to atmospheric effects. The mean values or the following 
parameters were computed from these measurements: (1) source SPL, (2) average Mach number, (3) liner impedance, 
(4) exit impedance, (5) static pressure, and (6) static temperature. Each of these mean values were used as inputs to 
the duct propagation codes, and the predicted attenuations were compared to the mean of the measured attenuations. 

Table 6 Comparison of measured and predicted liner attenuation for three selected codes using a segmented liner 


Mach 

kHz 

Measured 

C1I2DDS 

LEE2DDS 

CH3DPA 

0.00 

1,0 

9.52 

10.25 

10.43 

12,08 

0.00 

2.0 

7.44 

8.14 

8.23 

22,30 

0.00 

3.0 

2.03 

2.03 

2.07 

1.57 

0.25 

1.0 

8.55 

11.42 

11.25 

14.88 

0.25 

2.0 

23.54 

23.08 

24.08 

22.29 

0.25 

3.0 

6.06 

5.85 

5.74 

5.67 

0.40 

1.0 

5.57 

7.53 

3.81 

9.59 

0.40 

2.0 

14.44 

13.95 

13.93 

13.06 

0.40 

3.0 

7.95 

7.12 

6.99 

7.31 


Table 6 compares the predicted code attenuations (using the mean value f or the code input parameters) with the av- 
erage attenuation from the measurements. As observed in Table 6, the mean of the measured attenuation is in generally 
good agreement with code predictions with one exception. Note that the code that uses the parabolic approximation 
(i.e„ CH3DPA) does not agree well with the measurement at 2.0 kHz, and Mach zero. Among the most likely sources 
of the discrepancy between the CH3DPA result and measurement (i.c., at 2.0 kHz) are reflections from the leading or 
trailing edge of the liner, or the generation of w ide-angle modes in the duct. Both of these effects are neglected by the 
low fidelity code, CH3DPA. 
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B . E f ficien cy Assessment 


Two important metrics often used to determine llie efficiency of a computer code are wall clock lime mid computer 
memory requirements. The suite of codes presented in this study solve a discrete set of acoustic equations using cither 
a direct solve (c.g., CH2DDS, LEE2DDS, CH3DPA, CH3DDS, LEE3DDS) or an iterative solve (c.g., LEE2DIS, 
LEE3DIS) strategy to compute the acoustic solution. The codes that are based on an iterative solve strategy were 
designed for computers with low memory, whereas those that implement direct solve strategies are designed for high 
memory machines (with the exception of CH3DPA). Generally speaking, the wall clock time is not as important on 
low memory machines because they are often under the control of a single user and there is little competition for wall 
clock time from other users. On the other hand, liigh memory machines often support many users that compete for 
available wall clock time. In this section, we choose to place focus on how these codes are used at the NASA Langley 
Research Center (LaRC), where the most important efficiency metric for deciding on the choice of a code is the wall 
dock time. This occurs for three primary reasons: (1) eodes are often run on a massively- parallel computer within 
a multi-user environment where users must compete for wall clock time, (2) much emphasis is placed on impedance 
eduction, in which the selected propagation code must be run hundreds (sometimes thousands) of times to educe the 
unknown impedance of a test liner, and (3) InRC is working toward the ability to model full-scale, large commercial 
engine nacelles with high fidelity codes to get a quality prediction (high fidelity eodes require significantly more wall 
clock time than low fidelity eodes). Again, it is reemphasized that other research organizations may utilize these codes 
for other reasons or may have different computing facilities that may dictate choosing a different metric for analyzing 
code efficiency. 

Table 7. Wall clock time consumed on Columbia clusters by in-duct codes on baseline grid 


Code 

No. of 

Governing 

Numerical 

No. of 

Wall 

Name 

Dims 

Equations 

Method 

CPUs 

clock 

CH2DDS 

2 

Convec ted 1 lei mhol t z 

Cubic Finite Element 

8 

48s 

LEH2DDS 

2 

Linearized Euler 

Linear Finite Element 

8 

60s 

LEE2DIS 

2 

Linearized Euler 

Ps eudo Ti me- S tepping 

1 

210 s 

CH3DPA 

3 

Convected Helmholtz 

Parabolic Approximation 

1 

7s 

CH3DDS 

3 

Convected Helmholtz 

Cubic Finite Element 

32 

912s 

LEE3DIS 

3 

Lineari/ed Euler 

Ps eudo T i me- S teppi ng 

4 

5400s 


Table 7 gives the wall clock time on the baseline grid for the various eodes. The six columns in the table correspond 
to the code name, number of’ spatial dimensions, governing differential equations, the type of numerical method, the 
number of CPUs used, mid the wall clock time in seconds. Although the wall clock limes are those to obtain the 
solution in the hardwall duct, there was little change in the wall clock time when the wall was lined with acoustic 
material. As should be expected, the two codes that use direct solve strategies (CII2DDS, LEE2DDS) and run on more 
than one CPU require considerably less wall clock time than LEE2DIS, which employs an iterative solve strategy and 
is sequential. It is also noted that the codes using a direct solve strategy (CII2DDS, CII3DDS, LEE2DS) can be run in 
parallel to further reduce the wall dock time given in Table 7, although they do not scale linearly with the number of 
CPUs. Table 7 shows that the wall clock lime increases by a factor of 19 when the cubic finite element analysis (Le,, 
with a direct solve strategy) is increased from 2D (CH2DDS) to 3D (CH3DDS). Note that the number of CPUs had 
to be increased from eight to thirty -two to allow enough memory for the CH3DDS code to run (i.c., when compared 
to CH2DDS) mid a nineteen -fold increase in wall clock lime is still observed. In contrast. Table 7 shows that the 
wall clock time increases by a factor of 26 w hen the pseudo time- stepping code (i.c., with an iterative solve strategy) 
is increased from 2D (LEE2DIS) to 3D (LEE3DIS). It should be noted that increasing the number of CPUs for the 
LEE3DIS code beyond the four in fable 7 leads to an increase in wall dock time over that in Table 7. As expected, 
the lowest fidelity code in the suite, CH3DPA (CH3DPA; uses the parabolic approximation), is the most time efficient. 
The sequential code, CII3DPA, is observed to be more than two orders of magnitude faster than CII3DDS and nearly 
three orders of magnitude faster than LEF3DIS. 
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VII. Summary and Concluding Remarks 


A Suite of seven of the most commonly used duct propagation codes (three vortical and four non- vortical) has been 
evaluated in this paper. The suite consists of the most commonly used codes and methods. These codes support 
propagation through acoustically treated inlet and aft-fan ducts in the presence of mean flow and, collectively, represent 
the state -of-art in the liner industry. The following is a brief summary of w hat was achieved in this study: 

1 . The basic assumptions, governing differential equations, boundary conditions, and numerical methods underly- 
ing each code development have been quail tilled. 

2. Each code has been evaluated on the same problems, on the same computer, and using identical grids. 

3. The predicted attenuations in a hard wall duct were computed for nine combinations of frequencies and Mach 
numbers, and were compared to known values determined from analytical models, 

4. The predicted attenuations (for each code) over a section of an infinitely long duet with a uniform impedance 
liner along the Lop wall were computed for 36 combinations of Mach number, frequency, mid lining impedance, 
and were compared to values obtained from mode theory. 

5. As part of tills evaluation, attenuation and code input parameters were measured eight different times over a 
period of one year using the same conventional perforate- over ‘honey comb liner in the NASA Langley Grazing 
Incidence Tube (GIT). The measured code inputs were then averaged and used as input for code attenuation 
predictions* These predictions were compared to the measured mean attenuations for a 140 dB incident sound 
pressure level, for three flow rales, and source excitation f requencies of 1,0, 2,0, mid 3.0 kHz* 

The following primary findings have resulted: 


1. Wall clock time increases of an order of magnitude or more are observed for 3D codes relative to the corre- 
sponding 2D versions of the same codes. As should be expected, codes using mi iterative solve methodology are 
considerably more memory efficient and less time ef ficient than those based Upon a direct solve methodology. 

2. Attenuations predicted with each code are in good agreement with exact results for the hardwall duct. 

3. Attenuations predicted with the CH2DDS, LEE2DDS, CH3DDS, and CH3DPA codes over a portion of an 
in finitely -long duct with uniform acoustic treatment compare well with results from mode theory for the 36 
combinations of Mach number, frequency, mid wall impedance used in this study Of these Tour codes, the 
CH3DPA is considerably more time efficient. The LLE2DIS attenuation predictions compare well with mode 
theory at low Mach numbers, but tend to diverge from both mode theory and other code results as the flow 
Mach number is increased. Only the finite element simulation of the convected Helmholtz equation using cubic 
Ilermite polynomials w r as able to accurately predict extremely large liner attenuations. 

4. Attenuations predicted for a conventional perforate- over-honey comb liner installed in the GIT compare well 
with measured attenuations when compared on a statistical basis. 

This investigation has provided valuable insights with regards to strengths and weakness of several duct propagation 
codes and is expected to form the basis for continuing code development at NASA Langley Research Center, 
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